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We propose efficient measurement procedures for the self-energy and vertex function of the An- 
derson impurity model within the hybridization expansion continuous-time quantum Monte Carlo 
algorithm. The method is based on the measurement of higher-order correlation functions related 
to the quantities being sought through the equation of motion, a technique previously introduced 
in the NRG context. For the case of interactions of density-density type, the additional correlators 
can be obtained at essentially no additional computational cost. In combination with a recently 
introduced method for filtering the Monte Carlo noise using a representation in terms of orthogonal 
polynomials, we obtain data with unprecedented accuracy. This leads to an enhanced stability in 
analytical continuations of the self-energy or in two-particle based theories such as the dual fermion 
approach. As an illustration of the method we reexamine the previously reported spin-freezing and 
high-spin to low-spin transitions in a two-orbital model with density-density interactions. In both 
cases, the vertex function undergoes significant changes, which suggests significant corrections to 
the dynamical mean-field solutions in dual fermion calculations. 

PACS numbers: 71.10.-w,71.27.-|-a,71.30.-|-h 



I. INTRODUCTION 

Continuous-time quantum Monte Carlo solvers (for a 
recent review, see Ref. 1) have become an important tool 
for the study of the Anderson impurity model (AIM) and 
its multi-orbital generalizations, due to their accuracy, ef- 
ficiency and ability to treat arbitrary interaction terms. 
While the AIM plays a fundamental role in various areas 
of condensed matter physics, the continuous-time solvers 
have been developed and primarily applied in the con- 
text of dynamical mean- field theory (DMFT),^ which 
maps correlated lattice models to an appropriately de- 
fined quantum impurity model. 

Impurity models with multiple orbitals are important 
for two different reasons: (i) The combination of den- 
sity functional theory and DMFT'^ provides a formal- 
ism for the calculation and prediction of properties of 
strongly correlated materials. The description of materi- 
als with open d- or f-shells requires the solution of impu- 
rity models with up to five or seven correlated orbitals, 
respectively, (ii) In the context of cluster extensions of 
DMFT,^ the lattice is mapped onto a cluster of impuri- 
ties in order to account for spatial correlations, and it is 
desirable to solve clusters with as many sites as possible, 
in order to be able to infer reliable predictions for the 
infinite system. For cluster DMFT calculations of sim- 
ple models, the interaction expansion or weak-coupling 
algorithm, which is based on an expansion of the par- 
tition function in the interaction (henceforth referred to 
as CT-INT^) and the related continuous-time auxiliary 
field algorithm (CT-AUX^) have proven useful. Here 
the number of interaction terms and hence the perturba- 
tion order grow linearly with the number of cluster sites. 
In the context of ab initio calculations of correlated mate- 



rials, however, the number of interaction terms grows at 
least as the square of the number of orbitals and hence 
the hybridization expansion algorithm^^^^ (abbreviated 
CT-HYB) is the method of choice. 

For the latter, the calculation of the one-electron self- 
energy has proven problematic. The calculation from 
Dyson's equation, where it is evaluated as the difference 
between the inverses of two Green's functions, is highly 
susceptible to noise in the numerical data. In contrast to 
CT-INT, the Green's function in CT-HYB is not mea- 
sured as a correction to a known function (the noninter- 
acting Green's function Go), which leads to large noise 
in the intermediate to high-frequency region. Except for 
the low-frequency region, better statistics is needed for 
the calculation of the self-energy in CT-HYB to obtain 
results of comparable accuracy as in CT-INT. 

Similar problems arise in the calculation of the im- 
purity vertex function. While the vertex allows one to 
access response functions of a system, interest in this 
quantity has recently grown in particular due to the ad- 
vent of novel diagrammatic extensions of the dynamical 
mean- field theory. ^"^"^^ In these approaches, spatial corre- 
lations beyond DMFT are included through two-particle 
field theories, which involve the two-particle irreducible 
(in the dynamical vertex approximation^'*) or reducible 
(in the dual fermion approach*^'*^) vertex function. Thus 
far these approaches mainly relied on exact diagonaliza- 
tion (ED) or CT-INT for the calculation of the impu- 
rity vertex functions. Only few applications employing 
the CT-HYB algorithm for measuring the two-particle 
function have been reported because of the aforemen- 
tioned problems. The measurement of frequency depen- 
dent two-particle functions is particularly challenging for 
multiorbital models. For example, in a recent letter by 
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FIG. 1: (Color online) Typical segment configuration in 
the hybridization expansion continuous-time quantum Monte 
Carlo simulation of the one-orbital AIM. Each segment marks 
an imaginary-time interval in which an electron of given spin 
resides on the impurity. The overlap between the up- and 
down-spin segments (blue or gray shaded area) determines 
the interaction energy. 

Park et at, this problem was circumvented by neglect- 
ing the frequency dependence of the irreducible vertex 
function and hence setting it to a constant.^* 

On the other hand, the CT-HYB approach with its low 
perturbation order and favorable sign statistics^^ appears 
to be the most suitable, in principle, for measuring vertex 
functions in realistic multiorbital models. 

In this paper, we propose an efficient method for cal- 
culating the self-energy and vertex function which elim- 
inates the noise problems. Through the equation of mo- 
tion, the self-energy and vertex function are related to 
higher-order correlation functions, which can be mea- 
sured in CT-HYB. In combination with a recently devel- 
oped method to eliminate the Monte Carlo noise through 
a representation in terms of orthogonal polynomials,^^ 
the present approach yields considerably more accurate 
results compared to the standard measurements. 



II. METHOD 
A. Model 

We consider the single-site, multiorbital AIM with 
interactions of density-density type, described by the 
Hamiltonian 

^^AiM = ^yljk^ + ^^^^ U^jUinj 

ki i ij 

(1) 

where latin subscripts denote a 'flavor' index, i.e. a com- 
bined index for spin- and orbital degrees of freedom. The 
operators c| (cj) create (destroy) an electron with flavor 

i on the impurity site, while /j^^ (/ki) creates (destroys) 
a conduction band electron in band i with momentum k 
(ejj. is the corresponding dispersion) . The impurity levels 
are denoted by and the matrix Uij contains the vari- 
ous interaction parameters for the interactions of density- 
density type. Electrons with band-flavor index j and 



momentum k are allowed to couple to electrons with any 
other flavor index i as described by the hybridization ma- 
trix V^-' . Integrating out the noninteracting conduction 
band electrons gives rise to the hybridization function 

Aa.M=E^^^- (2) 

Here we use an implementation based on the segment 
picture of the hybridization expansion algorithm.^ The 
segment picture applies whenever operators of a given 
flavor appear in alternating order, which is always the 
case for an interaction of density-density type. 

In this case the trace over the sequence of impurity 
creation and annihilation operators (which defines the 
Monte Carlo configuration) does not have to be computed 
explicitly. In order to evaluate the ratio of traces one in- 
stead only needs to compute the length of the segments 
for the different flavors and their overlaps. This yields 
a very efficient algorithm for multiorbital problems with 
density-density interaction, which — as long as the deter- 
minant calculation dominates the computational effort — 
scales linearly in the number of flavors.^ A possible seg- 
ment configuration for the onc-orbital AIM is illustrated 
in Fig. 1. The segments represent imaginary-time inter- 
vals in which an electron of given flavor (spin) resides on 
the impurity. Overlapping segments correspond to time 
intervals in which the impurity is doubly occupied. As 
we will see below, the additional higher-order correla- 
tion functions which arise in the expressions for the self- 
energy and vertex function in the improved estimator can 
be measured at essentially no additional computational 
cost in this segment representation. The generalization 
to the case of interactions of non-density-density type, 
such as spin-flip and pair hopping terms, is in princi- 
ple straightforward. The measurement of the required 
correlation functions in the hybridization expansion al- 
gorithm, however, becomes more involved. 

In the following subsection, we first introduce and illus- 
trate the procedure for the self-energy and consider the 
generalization to the vertex function in a second separate 
subsection. 

B. Self-energy 

The asymptotic tail of the self-energy can be obtained 
from a high-frequency expansion, which requires mea- 
surements of single- and two-particle equal-time corre- 
lators (see e.g. Ref. 20 and references therein). This 
eliminates the noise at high frequencies. However, the 
choice of the cutoff value cannot easily be automatized, 
and as we will show (see, e.g., the inset of Fig. 6), the 
noise problem of the standard measurement is significant 
even at frequencies for which the self-energy clearly has 
not yet reached its asymptotic behavior. 

The method presented here yields accurate results 
over the entire frequency range and considerably re- 
duces the noise at intermediate to high frequencies. An 
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FIG. 2: (Color online) The 3! ways of connecting hybridiza- 
tion lines in a configuration with 3 segments for a timeline 
of a given flavor. The sign of each diagram is indicated. For 
non-diagonal hybridization, the lines also connect segments 
with different flavors. 



alternative approach to reduce the Monte Carlo noise 
is to measure observables in an orthogonal polynomial 
representation."'^^ This corresponds to an advantageous 
change of basis, which yields a compact representation 
of observables and allows one to filter the Monte Carlo 
noise, without any loss of information. We note that such 
a procedure does not reduce the required statistics and 
hence for optimal results the two approaches should be 
combined. 

The computation of the self-energy from higher-order 
correlation functions of the impurity model has been in- 
troduced previously in the context of numerical renor- 
malization group (NRG)2i calculations. Here we pro- 
vide the expression for the multiorbital case and demon- 
strate the usefulness of this scheme for the CT-HYB al- 
gorithm. An expression for the self-energy in terms of 
a two-particle correlator is obtained from the equation 
of motion for the Green's function. One may write the 
equation of motion in terms of a derivative with respect 
to either the first (t) or second (r') of its arguments. 
For reasons outlined in the appendix, we choose the sec- 
ond option. Here we only present a brief outline of the 
derivation. Details can be found in Appendix A. 

The time derivative of the Green's function 

G',fc(r-r'):--(T,c,(r)4(r')) (3) 

is given by 

dr'Gabir - r') = S{t - T')5ab - {TrCa{T)dr^cl{T')), (4) 

where is the usual time-ordering operator. Its equa- 
tion of motion follows from the one for the operator cj, 
taken in the Heisenberg representation: 



(5) 



We introduce the following correlation functions: 

GtAr-r'):=-{TrCa{T)fl{T')), (6) 
FLir r') := -{T^Ca{T)cl{T')n,{T')) , (7) 

together with their Fourier transforms G^^^f,(w) and 
F^j(w). With an application in the context of DMFT 
in mind, we furthermore introduce the noninteracting 
Green's function of the impurity model. 



where Aah(w) is the hybridization function. Evaluat- 
ing the commutator in Eq. (5) with the Hamiltonian (1) 
yields 



*jb 
k ' 



(9) 



SO that Eq. (4) in Fourier space can be written 
Gabiii^) = Go^abiiv) - ^G'Qi(w)Ajj(ii/)G'ojfc(w) 

(10) 



The function G'^^^{iv) in turn can be eliminated by ex- 
pressing it in terms of the impurity Green's function 
through its equation of motion. In Fourier space, it reads 



Gtab{^y)=Y.^a^{.^y)^ 



ib 



IV ~ el 



(11) 



Inserting this into (10) and using the expression for the 
hybridization function (2), the hybridization terms can- 
cel. Comparing the resulting expression with Dyson's 
equation (see Fig. 3), 

Gabiiv) = Go,Q&(ij^) + E Gai{w)'i^ij{iy)Go^jb{w), (12) 



we see that the impurity self-energy can be expressed in 
the following form: 

^ab(ll^) - J E Gah^^W.b + Ub,)Fi{lv), (13) 



which for a single orbital model reduces to the result 
given in Ref. 21. This equation relates the self-energy to 
two measurable quantities, the impurity Green's function 
Gab{ii^) and an additional correlation function F^f^{iv). 

In order to show how correlation functions are mea- 
sured in the CT-HYB, we recall that in this algorithm, 
one samples over configurations whose weight is given by 
a determinant of hybridization functions times a trace 
over a sequence of operators. In the segment picture, a 



h = a 



^Oafc(*^) = ~ £b)5ab - l^ab{w), 



(8) 



FIG. 3: Diagrammatical representation of the Green's func- 
tion Gab in terms the self energy E, Eq. (12). Thick lines 
denote fully dressed propagators, thin lines denote bare prop- 
agators. 
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FIG. 4: (Color online) Self-energy as a function of Mat- 
subara frequencies for the half-filled Hubbard model on the 
Bethe lattice (bandwidth 4£) as calculated from Dyson's equa- 
tion to illustrate the noise problem. The parameters are the 
same as in Ref. 12, U/t = 4 and T/t = 1/45. The results 
were obtained by measuring the Green's function directly on 
Matsubara frequencies and in imaginary time using 1500 bins 
and subsequent Fourier transform for the latter. The high- 
frequency tail lim^_+oo = U^{n){l — {n))/{iv) is shown 
for comparison. 
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FIG. 5: (Color online) CT-HYB data for the self-energy 
computed from Eq. (13) for the same model and parameters 
and measured within the same simulation as the results in 
Fig. 4. The Green's function and the two-particle correla- 
tor have been measured on the Matsubara axis. The line 
with closed symbols (blue) shows the result from an indepen- 
dent DMFT calculation using the interaction expansion algo- 
rithm (CT-INT), where the self energy has been obtained via 
Dyson's equation with the Green's function measured on the 
Matsubara axis. The noise problem at intermediate and high 
frequencies does not exist in the CT-INT and is resolved in 
CT-HYB using the improved estimator. The inset shows a 
blowup of the intermediate frequency range. 



configuration C is visualized by a collection of segments 
on the timeline from to /? for each flavor (/3 is the inverse 
temperature). A typical configuration for a single-orbital 
model with Hubbard interaction Un^n^ is depicted in 
Fig. 1. The start of a segment (open circles) is associ- 
ated with an impurity creation operator, while an impu- 
rity annihilation operator is associated with the segment 
endpoint (closed circles). A segment hence marks the 
imaginary time interval in which the impurity is occupied 
by an electron of a given flavor. The segments are con- 
nected by hybridization lines in all possible ways, which 
is the interpretation of the determinant (see Fig. 2). For 
the case of non-diagonal hybridization considered here, 
the segments are also connected by hybridization lines 
linking different flavors (not shown). 



A contribution to the Green's function of a particu- 
lar Monte Carlo configuration is obtained by cutting all 
hybridization lines connected to a given pair of a cre- 
ation and an annihilation operator, leaving a configura- 
tion with two unconnected operators. This corresponds 
to removing row i and column j from the matrix of hy- 
bridization functions A''. Denoting the resulting ma- 
trix as A!^j, the contribution of the particular configu- 
ration is essentially the ratio between the matrices af- 
ter and before removing the hybridization functions, i.e. 
(— 1)'+^ det Af^/det A*^. This ratio is precisely the 
element of the inverse of the matrix of hybridization func- 
tions, denoted by Mj^. Hence the Green's function de- 
fined on the interval from to /3 is measured according 



^9,10 

Gab{T~T') 



where 



J2 Mf,J-{T-T',T^~Tl)Sa.Jb,p) , (14) 
\al3=l I MC 



5-{t,t') := sgu{T')5{T - t' - e{-T')l3). (15) 



The only difference in the measurement for the function 
F^i^{t — t') is the presence of the additional density op- 
erator. Hence the measurement formula reads 



\a/3=l 



Tp)nj{Tl)5a,c,5b.fi 



MC 

(16) 



In the segment picture, the matrix element nj{Tp) (one or 
zero) of this operator is simply determined by examining 
whether or not a segment is present in the timeline for 
fiavor j at time rj . Hence this quantity can be measured 
at essentially no additional computational cost. Note 
that this function according to (13) only contributes if j 
is different from the index b (and /3) and therefore Uj is 
never evaluated at the position of the creator of flavor b 
at time t|. 

It is convenient to accumulate 

(GE)„,(r - r') = (1/2) Y,{U,b + UM{t - r') (17) 
j 
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FIG. 6: (Color online) Same as in Fig. 5, albeit with the 
Green's function and the two-particle correlator measured in 
the Legendre polynomial basis and transformed back to Mat- 
subara frequencies. This allows to efficiently filter the Monte 
Garlo noise. The inset compares the result at intermedi- 
ate frequencies to that of Fig. 4 obtained from the Fourier 
transformed imaginary time measurement (data are measured 
within the same simulation). 
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FIG. 7: (Color online) Comparison of Monte Carlo data for 
the imaginary part of the self-energy (obtained with the same 
method as in Fig. 6) with the exact result for a correlated 
site coupled to a single bath level. Some points of the Monte 
Carlo result are skipped to improve visualization. The result 
computed from Dyson's equation (with the Green's function 
measured on the Matsubara axis within the same simulation) 
illustrates the substantial improvement in accuracy. 



directly instead of the individual quantities F^^^It — t'). 
This is then analogous to the measurement of G times S 
in the CT-AUX algorithm.^ Note that these correlators 
can be measured in any basis by appropriately transform- 
ing the measurement rules, e.g. by taking the Fourier 
transform to measure directly the Matsubara coefficients. 

If the interaction is of non-density-density type, the 
segment picture has to be abandoned. In this case, the 
equation of motion generates additional terms, which re- 
quire the measurement of correlation functions of the 
form (rT-Ca(T)c[(r')cj(r')c/j(T')) (see appendix C). In ad- 
dition to computing the ratio of determinants, the mea- 
surement of such functions requires the evaluation of the 
ratio of traces over the operators in the particular con- 
figuration with and without the operators c^j{T')^C}^{T') 
inserted at the corresponding time. 

The major benefit of rewriting the self-energy in the 
form (13) instead of computing it from Dyson's equation 

T.{w) ^ Gaiiv)-^ - G{w)-^ (18) 

is that it is expressed in terms of a ratio of two mea- 
sured quantities instead of a difference of two functions. 
Forming the difference between an exactly known and an 
approximate quantity is susceptible to numerical errors, 
since, as already noted in Ref. 21, for a ratio only the 
relative error propagates, while in a difference, the ab- 
solute error (here the absolute error of G~'^) propagates. 
In DMFT, Go is computed from the self-energy of the 
previous iteration and is not known exactly. Error can- 
cellation however cannot be expected since Gq and G are 
computed in two independent Monte Carlo runs. Since 
G decays as l/{iv) it is clear that computing the self- 
energy according to Eq. (18) will lead to large errors in 
particular at intermediate to large frequencies. 

In the following we show results illustrating the ad- 
vantages of the improved measurement for the CT-HYB. 



We concentrate on the DMFT solution of the Hubbard 
model on the Bethe lattice with bandwidth 4i. The pa- 
rameters are the same as in Ref. 12. In order to ensure 
comparability, we have performed the various measure- 
ments within the same simulation, i.e. all quantities have 
been averaged over the identical sequence of Monte Carlo 
configurations (including the results of Fig. 4). 

Figure 4 shows the self-energy determined in the stan- 
dard way, i.e. by measuring the Green's function and 
calculating the self-energy using Dyson's equation. The 
Monte Carlo noise is clearly visible for intermediate to 
large frequencies, as anticipated. It is important to note 
that there is also considerable noise in the region where 
the self-energy has clearly not yet reached its asymptotic 
behavior (cf. also inset of Fig. 6 below). While the re- 
sults are similar for the frequency and imaginary time 
measurement at small to intermediate frequencies, the 
Monte Carlo error steadily grows with the frequency for 
the measurement in Matsubara frequencies. The noise 
sets in much below the Nyquist frequency (of approxi- 
mately Vn/t 50). In the imaginary-time measurement 
the noise is suppressed above the Nyquist frequency and 
the result smoothly approaches the tail, which is enforced 
in the calculation of the Fourier transform. 

Figure 5 clearly shows the advantage of the improved 
measurement: The error is considerably reduced over the 
entire frequency range. The number of measurements re- 
quired for a given accuracy at intermediate and large 
frequencies is hence greatly reduced. In the same figure 
we compare the improved measurement of CT-HYB with 
a result obtained from a comparable run of the CT-INT 
using Dyson's equation, with the Green's function also 
measured on the Matsubara axis. We note that this is 
meant as an illustration of the convergence properties at 
high frequency. Here we do not attempt the intricate task 
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of separating the various effects that influence the perfor- 
mance and accuracy of these complementary algorithms, 
which furthermore scale very differently with the number 
of flavors.^ Comparable refers to the fact that for both 
runs we have used similar simulation times that yield 
converged results in practice. A detailed performance 
comparison of the algorithms was presented in Ref. 12. 

The result illustrates that the noise problem does not 
exist in the CT-INT algorithm. The reason is that in CT- 
INT, the Green's function is measured as an 0{{l/ii')^) 
correction to the bare Green's function Gq. As a conse- 
quence, the evaluation of the self-energy is considerably 
more stable (see below). We find no advantage of calcu- 
lating the self-energy according to Eq. (13) in CT-INT, 
at least for the model considered here. 

We note that in principle one may evaluate the self- 
energy without explicitly measuring the Green's function 
G'ah(ii^). Combining Eqs. (12) and (13), one sees that the 
Green's function may be written in the form 



where the matrix Aan/ is defined as 

A , — A , 



(19) 



(20) 



The self-energy is obtained by multiplying (17) by the in- 
verse of (19). We find that this procedure yields results of 
similar but somewhat worse quality than the ones shown 
in Fig. 5, so that it does not reduce the effective compu- 
tation time for given accuracy. 

In Fig. 6 we show results obtained by combining the 
improved estimator with an efhcient method to suppress 
the Monte Carlo noise. The latter is based on a repre- 
sentation in terms of Legendre orthogonal polynomials. 
The transformation to the Legendre basis is exact, as 
in the Fourier case. Both G{iv) and F{iv) have been 
measured directly in the Legendre basis and transformed 
back exactly to Matsubara frequencies after the simu- 
lation. Appropriately choosing the cutoff in this basis 
allows to eliminate the Monte Carlo noise without loos- 
ing physical information (for details see Ref. 19). This 
method yields very accurate, noise-free results over the 
entire frequency range and captures the high-frequency 
tail correctly. We expect that Monte Carlo data mea- 
sured in this way will allow a more stable analytical con- 
tinuation. 

In order to show that the proposed method not only 
yields smooth, but indeed correct results, we compare 
the Monte Carlo data to an exact result. A suitable test 
case is a correlated impurity site coupled to a single bath 
level at e = for which the Hamiltonian (1) reduces to 



H = t[c^f + fc] + Un^n^ 



(21) 



const.. Although the hybridization function has no struc- 
ture in this case, this is not a trivial problem for the 
Monte Carlo approach, since all diagrams in principle 
have to be sampled to obtain the exact solution. The 
result at half-filling is shown in Fig. 7 for [//i = 4 and 
temperature T/t ~ 1/50. To improve the visualization, 
we have skipped some points from the calculated result. 
The curves lie on top of each other. That this is indeed 
not trivial can be seen from the result computed from 
Dyson's equation using a Green's function measurement 
on the Matsubara axis, which has again been accumu- 
lated within the same simulation. 



C. Vertex function 

To set the stage for the discussion of the vertex func- 
tion, we define the two-particle Green's function as 

Xabcd{Ta,Tb,Tc,Td) := {TrCa(Ta)cl{Th)Cc{Tc)cli{Ta)) . (22) 

Its Fourier transform is Xabcd{it^ai it^b, ii^c, ii^d)- With the 
disconnected part of the two-particle Green's function, 

Xabcdi'^^a, iVb, iVc, Wd) ■=GabiiVa)S^^,„^Gcd{Wc)S^^Ma 

-G ad{iv a)5y^,vaG cb{iv c)5y^,y^ 

(23) 

we obtain the connected part (cf. Fig. 8): 

XTbcd{i^a,il^b, iVc, ifd) -^XabcdiWa^Wb, Wc, W d) 

- X°abcd{i^a,Wb,iVc,Wd), 

(24) 

in terms of which the two-particle vertex function is de- 
fined as 

labcdij'^a,i>^b,iVc,Wd) := G^l'{i'^a)G^^,{Wc)^ 

X x''^b''c'd'{iVa,ivb,iVc,Wd)G^,l{wb)G^}^.i{wd). (25) 

The connected part of the two-particle Green's func- 
tion is the difference between a two-particle quantity and 
a product of two single-particle Green's functions. The 
calculation of this difference is susceptible to numeri- 
cal errors. Since single- and two-particle quantities have 
vastly different statistical errors, one cannot expect these 



d a 

+ 

c b 



d a ■ 



Xtrhcd 



Xtilcd 



and which is trivially solved by exact diagonalization. 
The hybridization function for this case is A(r) = t^/2 = 



FIG. 8; Diagrammatical representation of the two-particle 
Green's function Xabcd in terms of its connected part xTbik 
and its disconnected part Xabcd s-nd definition of the vertex 
function 7, Eqs. (23)-(25). The lines with arrows denote fully 
dressed propagators. 
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FIG. 9: (Color online) Real part of the spin-up-up compo- 
nent of the connected part of the two-particle Green's func- 
tion for a correlated impurity coupled to a single bath-level 
for T/t = 1/10 for fixed v' and bosonic frequencies uj = Q 
(circles), A-kT (triangles) and 2QnT (squares). The parame- 
ters are otherwise the same as in Fig. 7. Solid lines show 
the exact diagonalization data and open symbols the results 
from the improved Monte Carlo measurement (measured on 
the Matsubara axis). The usual Matsubara axis measurement 
based on Eqs. (22)-(25) is shown by crosses. The two results 
are hardly distinguishable on this scale. 



to cancel. The error in the vertex function itself is further 
amplified by multiplying the connected part with inverse 
Green's functions. This leads to deviations in particu- 
lar when G is small, i.e. for large frequencies or in an 
insulator, where G{ii>) extrapolates to zero for w — > 0. 

Hence we seek to express the connected part in a simi- 
lar fashion as was done for the self-energy using the equa- 
tion of motion. As shown in appendix B, the connected 
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FIG. 10: (Color online) Real part of the spin up-up com- 
ponent of the reducible vertex function "i^^i^ computed from 
the connected part shown in Fig. 9. The standard compu- 
tation (Eqs. (22)-(25), crosses connected by dashed lines) 
is clearly numerically unstable for large values of v^, while 
the improved estimators yield the proper asymptotic behav- 
ior. The improved estimator computed from Eqs. (25), (26) 
(closed symbols) and from Eqs. (25), (31) (open symbols) 
yield similar results. They also yield more accurate results 
even before the vertex approaches its asymptotic behavior and 
for small frequencies Vn (clearly visible for cu = 20nT). 



part of the two-particle Green's function can be written 
in the alternative form 

- ]^^^{Uji + Uij) Fl^{iVa) Xtbcd{iVa,iVb,iVc:iVd) 

(26) 

A priori it is not clear that computing the connected part 
according to this expression has any advantage compared 
to Eq. (24), since it still has the form of a difference of 
products of correlation functions. However we find that 
it indeed yields more accurate results. 

In addition to the single- and two-particle Green's 
function this expression involves the correlation functions 
^afc(*^) ^'^'-^ Fourier transform of 

Hii,cdi'ra,n,Tc,Td) {Trnj{Ta)Ca{Ta)cl{n)CciTc)cl{Td)). 

(27) 

In order to see how this function is measured, recall that 
in CT-HYB the two-particle Green's function with its 
arguments varying in the interval from to /3 is measured 
according to^ 

\abcd(.'^ab: ^cd: '^ad) 
/ '^'^ 



X S~{Tcd,T^ - Tg)5^{Tad,T^ " Tg)6a,aSb,l3Sc,',Sd,S 



MC 

(28) 



where is defined as before (15) and (S+(r, r') :— 
6{t — r' — 9{—t')P). Due to time translation invariance, 
the two-particle Green's function needs to be measured 
as a function of three independent time differences only. 
These have been chosen such that x is antiperiodic in 
Tab and Tcd^ while it is periodic in Tad- When taking the 
Fourier transform, the time differences Tab and Ted are 
associated with fermionic Matsubara frequencies u, v' ^ 
while Tad is associated with a bosonic frequency to. Note 
that the relation to the representation with four frequen- 
cies is such that 

Xahcd(y, v' , uj) EE Xabcd{w + lUJ , iv , w', iv' + iw). (29) 

The measurement formula for the correlator H^abcd r^ads 

^lbcdi^ab,Tcd,Tad) = 



U - MLM^p^)n, {t:)S- {Tab, <-Tf,)x 



X 5 {Tad, t'L - Tg)d+{Tad, " Tg)Sa,aSb,l3Scn^d.,S 



MC 

(30) 
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FIG. 11: (Color online) Spin-up-down component olthe real 
part of the reducible vertex function for a correlated impurity 
coupled to a single bath-level. The parameters are the same 
as in Fig. 7. Lines show the exact diagonaHzation result and 
open symbols show the improved Monte Carlo measurement. 
Black symbols correspond to the usual Matsubara axis mea- 
surement based on Eqs. (22)-(25) (note that the measurement 
in the Legendre basis leads to systematic deviations rather 
than the irregular noise observed in the Matsubara measure- 
ment). 



in analogy to F^^^. The correlator H''^^^^ can thus be mea- 
sured at essentially no additional computational cost to- 
gether with the two-particle Green's function. 

Single-particle quantities can be measured directly in 
imaginary time with essentially arbitrarily fine resolu- 
tion, since the number of imaginary time bins does not 
influence the performance of the algorithm. For two- 
particle quantities, which depend on three independent 
time differences, this is not the case, due to memory re- 
strictions. Hence we measure these quantities directly in 
the Matsubara or the Legendre basis. Instead of mea- 
suring the correlator {v, i^' ,uj) for all flavors j, we ac- 
cumulate {l/'2)J2j{Uja + Uaj)Hl^^^{iy,iy',u!) in order to 
save memory. However, this still requires one to measure, 
in addition to the two-particle Green's function, an ob- 
ject which is of the same size. In order to avoid this, one 
can replace Xabcd on the right hand side of Eq. (26) by 
the sum of its connected and disconnected parts, which 
leads to 



1 



a' I 



(31) 

where Aaa' is defined as in (20). Note that here the 
disconnected part instead of the full two-particle Green's 
function appears on the right hand side. In the following 
we assess the quality of the results obtained from the 



different measurement formulas for the vertex function. 

As a test, we again consider a correlated impurity hy- 
bridized to a single bath level. This allows us to compare 
Monte Carlo data to exact results. Figure 9 shows the 
real part of the spin-up-up component of the connected 
part of the two-particle Green's function, Rex]^^™Ji', as 
a function of the fermionic Matsubara frequency Un for 
fixed u' = 9nT {T/t = 10, U/t = 4). Circles, squares and 
triangles correspond to the bosonic frequencies w = 0, 
AttT, and 2QttT, respectively. The improved Monte Carlo 
measurements (open and closed symbols), as well as the 
usual Matsubara-axis measurements (crosses) agree very 
well with the exact diagonalization result (lines). The 
problems caused by the stochastic noise are not evident 
at the level of the two-particle Green's function, but only 
at the level of the vertex. 

The real part of the spin-up-up component of the 
reducible vertex is shown in Fig. 10. Here, due to 
the multiplication with inverse Green's functions, the 
noise in the data obtained from the standard Matsub- 
ara measurement grows considerably at large frequen- 
cies, while the improved measurement reproduces the 
correct high-frequency behavior. Both improved estima- 
tors, Eqs. (25), (26) and Eqs. (25), (31), yield results of 
comparable accuracy. In the case of larger bosonic fre- 
quencies, where the connected part becomes small, the 
improved accuracy of the new measurement procedure is 
evident even at the lowest Vn- 

We found that using the vertex function obtained via 
the improved estimators enhances the stability of dual 
fermion calculations. This is particularly the case for 
calculations in the insulating phase, where the improved 
estimators appear to be significantly more accurate and 
where it is otherwise difficult to obtain sufficient statistics 
due to the low perturbation order. 

Figure 11 shows the real part of the spin- up-down com- 
ponent of the reducible vertex at low temperature and for 
the same parameters as in Fig. 7 {T/t = 1/45, U/t = 4), 
measured in the Legendre basis. Again, one sees that 
at large bosonic ujm, deviations between the exact result 
(lines) and the standard Matsubara measurement (black 
symbols) appear already at small Vm while the improved 
estimators (symbols) yield more accurate data. 



III. APPLICATION 

A. Two-orbital model 

As an application of the methods described above, 
we calculate the self-energy and vertex function for a 
two-orbital model within the single-site DMFT. We con- 
sider the Hubbard model on the Bethe lattice with semi- 
elliptical density of states with bandwidth At. First, 
we re-examine the spin-freezing transition reported in 
Ref. 22 for a three-orbital model, and show that qualita- 
tively the same physics is found in the two-orbital model, 
with the interaction restricted to density-density terms. 
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FIG. 12: (Color online) Phasediagram in the space ol density 
(n) and interaction U ol the two-orbital Hubbard model on 
the Bethe lattice with J/U = 1/6 and T/t = 0.02. The thick 
black lines mark regions ol Mott insulating behavior. Circles, 
squares and stars, respectively, mark the parameters lor which 
the calculations in Figs. 15, 16 and 17 have been performed. 

We will work with the Hamiltonian (1), with the inter- 
action part explicitly given by 

+ {W - J)Y,ni,an2,a, (32) 

(J 

where a is the orbital index, a denotes spin, U and U' 
are the intra- and inter-orbital Coulomb interaction pa- 
rameters, J is the Hund's rule coupling coefficient and 
U' = U- 2 J. 



B. Spin-freezing transition 

The phasediagram in the plane of filling (n) and in- 
teraction U is shown in Fig. 12 for J/U — 1/6 and 
temperature T/t = 0.02. It reproduces the qualitative 
features reported in Ref. 22 for the three-orbital model: 
For small values of the density and interaction, we find 
a Fermi liquid phase, while for larger values, the model 
exhibits a frozen moment phase. The frozen moment 
phase is characterized by a spin-spin correlation function 
Css{t) ■= {Sz{t)Sz{0)) that approaches a non-zero con- 
stant at large times, as shown in Fig. 13. In this phase, 
the spin-spin correlation function at t — 1/(2T) (which 
we refer to as C1/2), is expected to be independent of tem- 
perature. In a Fermi liquid at low temperature, on the 
other hand, the spin-spin correlation function behaves as 
Css{t) ~ {T/ sin(7rTr))^ for times t sufficiently far from 
r = or 1/T, respectively, so that C1/2 ~ T^. In Fig. 14 
we show the ratio Ci/2(T = 0.02t)/Ci/2(T = O.Olt), as a 
function of filling, which clearly confirms this behavior. 
The ratio changes from the value 4 expected in the Fermi 
liquid phase to 1 expected for frozen moments, passing 
through 2 near the spin-freezing transition, indicating a 
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FIG. 13: (Color online) Imaginary-time dependence of 
the spin-spin correlation function for U/t = 8 and densities 
(n) = 1.90,1.72,1.53,1.29,1.10 and 0.72 (from top to bot- 
tom). For small densities, in the Fermi liquid regime, the 
spin-spin correlation at r = l/(2r) is proportional to T and 
hence approaches zero as T — >■ 0, while for high densities the 
spin correlation function approaches a constant, indicating 
the presence of frozen moments. 



T-linear behavior. Hence, in the vicinity of the transi- 
tion, the spin-spin correlation function decays unusually 
slowly, CssiT) ^ ^/^i consistent with the findings in 
Ref. 22. 

We note that a similar phasediagram for a two-orbital 
model, which plots the "quasi-particle weight" Z in the 
space of filling and interaction strength, has recently been 
reported in Ref. 23. Our spin- freezing transition line 
seems to resemble the contour-lines for fixed Z in Ref. 23, 
although at the temperature T/t — 0.02 of our calcula- 
tion, the strong deviations from Fermi-liquid behavior 
near the transition mean that in this regime a quasi- 
particle weight cannot be properly defined. It is an in- 
teresting open question whether and how the properties 
of a (strongly renormalized) Fermi-liquid are recovered 
as T ^ 0. 

The development of frozen moments is accompanied 
by a simultaneous change in the low-frequency behavior 
of the self-energy. In Fig. 15 we plot the imaginary part 
of the self-energy as a function of Matsubara frequen- 
cies for U/t = 8 and the densities indicated by circles 
in the phasediagram in Fig. 12. For small densities, the 
imaginary part of the self-energy ImS(w„ — >■ 0) extrap- 
olates to zero as expected for a Fermi liquid. In the 
presence of frozen moments, however, the electrons are 
expected to be scattered by these moments, resulting in 
a non-zero value — ImI](ii/„ — > 0) = rsgn(t'„). This 
behavior is evident from the imaginary time data. We 
find that in the transition region the low-frequency be- 
havior of the self-energy is well described by a power-law 
ImI](i/„)/i — j{i'n/t)"- The phase boundary in Fig. 12 
corresponds to a "square-root" frequency dependence, 
i. e. a = 0.5. Note that the horizontal axis is (vn/t)'^'^, 
so that a = 0.5 corresponds to a linear increase at low 
frequencies (black line in Fig. 15). 
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FIG. 14: (Color online) Ratio of the spin-spin correlation 
function Css{t = l/(2r)) for temperatures T = 0.02t and 
O.Olt as a function of density. Vertical lines mark the ap- 
proximate location of the transition, where the self-energy is 
consistent with a (vn/t)" behavior with a ~ 0.5. On the 
Fermi liquid side, obviously Css{r = l/(2r)) ~ T^, while 
Css{t ~ l/(2r)) — const, in the frozen moment phase. 



It is instructive to compute the self-energy on the real 
axis. Analytical continuation is a delicate issue, in par- 
ticular for Monte Carlo data, due to the presence of sta- 
tistical noise. Given the high quality of the Matsubara- 
axis data of Fig. 15, which have been obtained using the 
combination of improved estimator and Legendre noise 
filter, we nevertheless attempt to perform an analytical 
continuation using Fade approximants.^'* The results are 
shown by thick red (gray) lines in Fig. 16. The three 
values of the density correspond to the positions marked 
by squares in the phasediagram in Fig. 12. We first ob- 
serve that for (n) = 0.72, the real axis self-energy has 
the expected Fermi liquid form and is well described by 
the expression S(w) ^ k + Xui — i^uj'^ over a wide energy 
range (thin solid line) . For comparison we have also plot- 
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FIG. 15: (Color online) Imaginary part of the self-energy for 
various fillings for the two-orbital model at U/t = 8, obtained 
using the improved estimator measured in the Legendre basis. 
The positions of the corresponding parameters are marked by 
circles in the phasediagram in Fig. 12. The solid line is a fit 
proportional to (vn/t)", with a = 0.49. 
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FIG. 16: (Color online) Self-energy on the real axis obtained 
from the Matsubara axis data by analytical continuation using 
Fade approximants. Thick lines (red, gray in print) show 
the result obtained from the improved estimator measured 
in the Legendre basis. The dashed and dotted lines show 
results obtained from the improved estimator and Dyson's 
equation, respectively, measured in Matsubara basis. The 
positions of the parameters in the phasediagram are marked 
by squares in Fig. 12. In the Fermi liquid phase (left panel) 
the data obtained from the improved estimators measured in 
the Legendre basis is well described by the Fermi liquid form 
T,{uj) ~ {l — l/Z)Lo — iiiuj'^ over a wide energy range (thin solid 
lines), while the data obtained from the two other methods 
exhibits spurious features. 



ted the results obtained from the self-energy measured in 
Matsubara representation, using Dyson's equation (dot- 
ted lines) and the improved estimator (dashed). In par- 
ticular the self-energy obtained from Dyson's equation 
shows spurious features and a departure from the ex- 
pected Fermi liquid behavior already for energies close 
to the Fermi level. 

By construction, the Fade analytical continuation 
yields the most accurate result for small u. As a con- 
sistency check, we performed a least-squares fit of the 
Fade data to the Fermi liquid form in the energy win- 
dow [—0.1 : 0.1] (thin solid line), finding fi = 1.24 and 
A = (l — l/Z) — —1.0899, which corresponds to a weakly 
renormalized Fermi liquid with Z = 0.479. The value of 
A agrees remarkably well with the value obtained from a 
second-order polynomial extrapolation of the self-energy, 
A = —1.0845 for the improved estimator measured in 
Legendre (A = —1.0904 for the self-energy obtained from 
Dyson's equation). The deviation is less than 0.5%, while 
the fit of the Fade analytical continuation of the self- 
energy obtained from Dyson's equation already disagrees 
by more than 7%. Wc note that the Monte Carlo error 
on the self-energy is of the order of 2 • 10""*. A sum- 
mary of the parameters extracted from the Fade analyti- 
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Leg. + impr. est. 
AS #df X a \ 



impr. est. 



Dyson 



A 



0.02 -1.09052 - -1.08322 - -1.20196 - 

0.1 12 -1.08994 6.0-10"^ -1.08667 8.6-10"^ -1.16301 1.9-10"'' 
0.25 32 —1.0857 ^ 7.in — _i im /l.in^^ i nc/i n q in-3 



0.5 64 -1.070 3.5-10 



3.7-10"'' -1.103 

-3 



-1.11 



5.4-10" 
1.4-10" 



-1.064 9.3-10" 



-1.03 



1.2-10 



-2 



TABLE I: Coefficient A — 1 — 1/Z obtained from a feast- 
squares fit of ttie Pade data in the energy window [—AE : AE] 
using #df degrees of freedom. Tfie vaiues shouid be com- 
pared witti ttie result obtained from a second-order polyno- 
mial extrapolation of the imaginary part of the Matsubara 
self-energy, which yields A — —1.08449. 



els) and inter-orbital (right panels) components. The 
vertex is essentially featureless on the Fermi liquid side 
((n) — 0.72, red circles), but develops structure as the 
spin-freezing transition is approached ((n) = 1.10, green 
triangles), and notably in the the frozen moment phase 
{{n) ^ 1.72, blue squares). In particular the spin-channel 
(closed symbols) is strongly enhanced in the frozen mo- 
ment phase. 

One can anticipate that these features will induce sig- 
nificant changes when the vertex is used to calculate the 
momentum dependent self-energy in diagrammatic ex- 
tensions of dynamical mean-field theory. 



cal continuation for the various measurement procedures 
is given in Table I. 

The self-energy for (rt) = 1.13 (middle panel) shows 
incipient non-Fermi liquid behavior. Both the real and 
imaginary parts develop a peak around the Fermi level, 
qualitatively consistent with the appearance of a "square- 
root" non-analyticity near the spin-freezing transition. 
While the results from the improved estimators agree 
fairly well, the result obtained using Dyson's equation 
deviates rather strongly. Finally, for (n) = 1.53 the self- 
energy clearly exhibits non-Fermi liquid behavior, with 
even more pronounced peaks around the Fermi level. The 
non-zero value of the imaginary part at zero frequency is 
again consistent with the value obtained from a polyno- 
mial extrapolation of the Matsubara self-energy. As a fi- 
nal illustration for this model, we calculate the reducible 
vertex function according to Eqs. (25)-(26). Figure 17 
illustrates the evolution of the vertex for fixed i^' ~ — ttT 
and bosonic Matsubara frequency oj = 2ttT across the 
spin freezing transition, both for intra-orbital (left pan- 
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FIG. 17: (Color online) Real and imaginary parts of 
the reducible vertex function 'y^^'^ for the two-orbital Hub- 
bard model for v' = —nT and u = 2-kT across the spin- 
freezing transition. The left panels shows intra- and the 
right panels interorbital components of the spin (closed sym- 
bols) and charge (open symbols) channels. Circles (red), 
triangles (green) and squares (blue) correspond to densities 
(n) = 0.72, 1.10 and 1.72, respectively. The position of these 
parameters in the phasediagram is marked by stars in Fig. 12 
(in the corresponding color). The y-axis range in the upper 
left panel has been restricted to improve visualization. 



C. High-spin to low-spin transition 

As a second application, we study the same model at 
half-filling and additionally consider a crystal-field split- 
ting term A^f, defined as 



"-2,, 



(33) 



in Eq. (1). This model has recently been considered 
as a minimal model for the physics of LaCoOs.^^ The 
condition for half-filling is found by taking the chem- 
ical potential /i to be half the sum over the elements 
of any row or column of the interaction matrix, which 
gives ^1/2 = 3J7/2 — 5 J/2. For this model, a high-spin 
to low-spin, and associated orbital polarization transi- 
tion have been reported. To quantify the effect of our 
density-density approximation, we compute the phasedi- 
agram in the space of crystal field splitting and interac- 
tion strength. Comparison of the result shown in Fig. 18 
with the phasediagram for the rotationally invariant in- 
teraction (including spin-fiip and pair-hopping terms) re- 
ported in Ref. 26 shows that the differences are rather 
marginal and that all qualitative features are reproduced 




FIG. 18: (Color online) Phasediagram in the space of crystal 
field splitting Acf and interaction strength U of the half-filled 
two-orbital Hubbard model on the Bethe lattice with band- 
width At, J/U =1/4 and T/t = 0.02. The parameters for 
which the vertex function in Fig. 20 has been calculated are 
marked by stars. 
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FIG. 19: (Color online) Total moment (S'^), moment of or- 
bital 1, (5*1^) and density of orbital 1, (ni) as a funciton of 
U for various values of the crystal field splitting Acf. From 
left to right (and top to bottom) the curves correspond to the 
values Acf = 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5. 



by the density-density approximation. In particular we 
find a high-spin Mott insulating phase at large interac- 
tion and small crystal field splitting, and a low-spin or- 
bitally polarized insulator at small interaction and large 
crystal field splitting. These two qualitatively distinct 
insulating phases are separated by a metallic phase with 
an end-point near Acf/t = 4.5. 

The nature of the two insulating phases becomes ev- 
ident if one plots the spin expectation value (total mo- 
ment {SD or moment of orbital 1 {Sf^)) and the oc- 
cupancy of orbital 1 {{ni)) as a function of interaction 
strength. Several such traces corresponding to fixed val- 
ues Acf = 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5 are shown in Fig. 19. 
In the orbitally polarized insulator, the occupancy of or- 
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FIG. 20: (Color online) Real and imaginary parts of the 
reducible vertex function 'y^^^i for the two-orbital Hubbard 
model for v' = — vrT and uj = 2-kT at ^cf/t = 3.0 for the 
points indicated by stars (in corresponding color) in the phase 
diagram of the low-spin to high-spin transition in Fig. 18. 
The left panel shows intraorbital (for orbital 2) and the right 
panel interorbital components of the spin (closed symbols) 
and charge (open symbols) channels. Circles (red), trian- 
gles (green) and squares (blue) correspond to U/t = 6.0,7.0 
and 8.0. The corresponding points in the phasediagram are 
marked by stars in Fig. 18 (in corresponding color). 



bital 1 (which is shifted up in energy) is very low. The 
filling continuously increases with U through the metallic 
phase and eventually jumps to (rii) « 1 at the transition 
into the high-spin Mott insulating phase. A similar be- 
havior is seen in the traces for {S'^). 

We finally show in Fig. 20 the evolution of the reducible 
vertex function in the metallic phase, as one moves from 
the phase boundary to the orbitally polarized insulator 
to the phase boundary with the high-spin Mott insula- 
tor. The three data points, corresponding to Acf/t = 3.0 
and U/t = 6, 7 and 8 are marked in the phasediagram 
(Fig. 18) by star symbols. While rather featureless near 
the phase boundary to the orbitally polarized insulator, 
the vertex develops structure as the moment (S'f) in- 
creases (cf. the arrows in the upper panel of Fig. 19). 



IV. CONCLUSIONS AND OUTLOOK 

We presented efficient measurement procedures for the 
self-energy and vertex function within the hybridiza- 
tion expansion CTQMC approach. The improved esti- 
mators are based on higher-order correlation functions, 
which can be particularly easily measured in single- 
site (multi-orbital) models with density-density inter- 
actions, but also for generic models. In combination 
with a recently proposed noise-filtering scheme (Legen- 
dre representation), we were able to completely elimi- 
nate the noise problems at intermediate to high frequen- 
cies, which have plagued the "standard" evaluation of 
the self-energy and vertex function using the hybridiza- 
tion expansion method. With the noise-problem solved, 
one can fully exploit the performance advantages of CT- 
HYB in the strong-correlation regime and in the case of 
(single-site) multiorbital models. The accurate measure- 
ment of two-particle Green's functions and vertex func- 
tions should enable dual fermion^^'^^ or dynamical vertex 
approximation^^ calculations of multi-band systems, and 
thus (in combination with band-structure input) the ab- 
initio simulation of transition metal compounds which 
capture the effect of non-local correlations. While we 
have presented results for two-orbital models, the com- 
putational effort for the measurement of the vertex func- 
tion scales as the square of the number of orbitals (for 
diagonal hybridization) and thus simulations are feasible 
even for five-orbital models on a small computer cluster. 

We have further demonstrated the efficiency of the im- 
proved measurements with self-energy and vertex data 
for a two-orbital model with density-density interactions. 
In particular, we have revisited the phenomenon of the 
"spin freezing" transition, which was originally discov- 
ered in a three-orbital model, but which manifests itself 
in a very analogous manner also in the two-orbital case. 
In combination with the results by Ishida and Liebsch^^ 
for five-orbital models, this establishes the spin-freezing 
phenomenon as a generic feature of multi-orbital models 
with large Hund's rule coupling. This phenomenon leads 
to strong non- Fermi liquid effects in certain (U, n) regions 
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of the paramagnetic metallic phase. 
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Appendix A: Self-energy 

Following Ref. 21, we derive the expression for the 
self-energy for the Hamiltonian (1). The single-particle 
Green's function is defined as 

Gab{T - r') := -(r,c,(T)4(r')). (Al) 

In order to arrive at the equation of motion, for reasons 
outlined in appendix B, we take the derivative with re- 
spect to its second argument, 

b(T ~ r')(c,(r)4(r')) + <5(r - T'){S,ij')c,(T)) 
- Bir - r'){ca{r)dr4(r')) + 9{t' - r)(9.,4(r')c„(r)> 
= <5(r- r')(K(T),4(r')}) - (T.c„(r)9.,4(r')> 

= 6{t - T')6ab - {TrCa{T)dr-cl{T')). (A2) 

The equation of motion for the creation operator is 

dr,cl{r') = [H,cl]{r') (A3) 

and the commutator [H,cl] with the Hamiltonian (1) 
evaluates to 

i ij kz^' 

= ^bcl + \ J2{U,, + U,,)cln, + Y: (A4) 

Note that the commutator in (A4) involving the interac- 
tion term generates two terms with different order of the 
operators cl and rij. These can be commuted since the 
interaction matrix for identical flavors is zero: Ujj = 0. 
Inserting this expression into (A2), we obtain 

dr'Gabir - t') =S{t - T')5ab + G ab{T - T')eb 

+ \T.^U,, + U,,)Fi,{T-T') 
j 

+ Y.GtA^-r')V:'\ (A5) 

kj 



where we have used the above definition of Green's func- 
tion and introduced the functions 

Gtab{r-T') ■.^-{TrC,{T)fl{T')) 

Kbir-r') ■.^-{TrCa{r)cl(r')n,{T')). (A6) 
Using the definition of the Fourier transform 

G{iiy)^- dr dr' G{t - T')e"'^^-^ \ (A7) 
P Jo Jo 

(A5) can be written in Fourier space as 

GabMit,. Sb) =Sab + E Gt.i^l^W:^' 

kj 

+ I T.^U,b + Ub,)Fi,{iv). (A8) 
i 

Note that here no summation over repeated indices is im- 
plied (unless explicitly indicated). In DMFT it is custom- 
ary to define the bare Green's function of the impurity 
model as 

Galb(^^) = (^'^ - £b)5ab " ^ab{iv), (A9) 

where Aafc(w) is the hybridization function. Subtracting 
'YliiGai{iv)^ib{iv) on both sides of (A8), we can rewrite 
its left hand side as 

Y ^aiiii^) [{ii' - £b)Sib - Ajf,(ii^)] = Y Gai{i'^)GQ}b- 

i i 

(AlO) 

Hence multiplying both sides of (A8) by Go from the 
right using matrix multiplication, the equation of motion 
becomes 

Gab{w) = Go^ab - Y '^ai{iv)/^ij{iv)GQjb{w) 

+ EGka.MK''G0,.(*^) 
kij 

ij 

(All) 

The Green's function Gj^^j(r — r') in turn is determined 
from its equation of motion, 

dr.Gtbir r') = -{TrC,{r)dr, fUr')) , (A12) 
where now 

dr'fl{r') = [Hjl]{r'). (A13) 
The commutator is 

[H. flb\ = E ^k'K'. fL] + E K^'[4/k'„ /^j 

= <flb + Y4v^- (AM) 
i 
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Inserting this back into (A12) yields 

dr>Gt,{T r') = Gt,{r - r')4 + E ^a^ir - r')K'- 

(A15) 

In Fourier space, this becomes 

Gtbi^^) = E Ga.{^v)^^^. (A16) 

As can be shown by integrating out the bath fermions 
from the Hamiltonian (1), the hybridization function for 
this model reads 

Hence by inserting (A16) into (All) on sees that the 
terms involving the hybridization function cancel exactly: 

(A18) 

Comparing this with Dyson's equation, 

Gab{iv) = G^ahiw) + E Gaiiivy^.ijiiv^Gajhiiv), 

(A19) 

we can finally identify the self-energy as 

SabM = ^ E G-'i^'^mb + U,i)F^{tv). (A20) 



Appendix B: Vertex function 

In order to derive an analogous expression for the con- 
nected part of the two-particle Green's function or the 
impurity vertex function, respectively, we make use of 
the following operator identity (see e.g. Ref. 29): 

drATrA{Ta)B{n)G{T,)D{Td)) = 
{TrdrAija)B{n)C{T,)D{Td)) 
+5{Ta - n){Tr[A{Ta).B{Ta)]±G{T,)D{Td)) 
+6{Ta - T,){Tr[A{Ta),C{Ta)]±D{Ta)B{n)) 
+6{Ta - Td){Tr[A{Ta),D{Ta)]±B{n)C{T,)), (Bl) 

where is the (anti-)commutator [A, :— ABzL 

BA depending on whether one deals with fermionic or 
bosonic operators. For the two-particle Green's function 
defined as 

Xabcd{Ta,n,Tc,Td) := {TrCa{Ta)cl{n)Cc{Tc)cl{Td)) , 

(B2) 

we have 

dr^XabcdiTa, U, Tc, Td) =9r„ {TrCa{Ta)cl{n)Cc{Tc)c\{Td)) 
= {TTdr^Ca(Ta)cl{n)Cc{Tc)c\{Td)) 
- S{Ta - n)Gcd{Tc - Td)5ab 
+ ^(ja - Td)Gcb{Tc - n)6ad- 

(B3) 

The last two terms stem from the discontinuities of 
Xabcd{Ta,Tb,Tc,Td) at r<j = Tfc and Ta ^ Td. Defining the 
Fourier transform 



I rl3 rP rP 

T[f{Ta,n,Tc,Td)]:= - dTa du dr^ dTdf{Ta,n,Tc,Td)e"'''^'^e-"'''^''e"'''^'e-"'-'^'' = f{ii'a,iiyb,ii^c,ii^d), 
P Jo Jo Jo Jo 

(B4) 

I 



one has 

J^i^iTa - Tb)Gcd{Tc - Td)] = (iSy^ ,UbG cd{iv c)5u, ,Ud^ 
^[Xabcd{Ta, Tb, T^, Td)] = Xabcd{il^a, Wb, ^c, W d) ■ 

(B5) 

Using the equation of motion for Cairo), i.e. dr^Ca(Ta) — 
[H,Ca]{Ta), we nccd the commutator 

[H, Ca] = -EaCa - ^ E^^^" + Uaj)njCa - E K'^/kj' 

(B6) 



which leads to 

{iVa - £a)Xabcdiil^a, i^b, Wc, W d) = 

P {Gcd{iVc)&u^,vt^u,Mi5ab - Gcb{T'Vc)5iu^,tvJiu^^,ii^Jad) 

+ E ^k'' ^ijbcdiiva,ivb, ivc, ivd) 

kj 

+ 2 E^^J'^ + Uaj)Hlbcd(^'^a,Wb,Wc,il'd)- (B7) 

j 

Again, no summation over repeated indices is implied 
unless explicitly indicated. Here we have introduced the 
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correlation functions 

Xka6cd(^a,T6,Tc,rd) := (T^/ka(ra)4(T6)Cc(Tc)4(Td)), 

Hibcdi'^a,n,Tc,Td) ■■= {Trnj{Ta)Ca{Ta)cl{n)CciTc)c]^iTd)). 

(B8) 

The first one can again be eliminated by considering its 
equation of motion, 

^Ta Xkabcd iTa,Tb,Tc,Td) = 

{TrdrJUTa)cl{n)Cc{Tc)c\{Td)). (B9) 

Using the equation of motion for f^aiTa), 

drJUra) = -eUUra) - K^'cjira), (BIO) 

j 

one finds 

dfccc / \ Of fccc / \ 

- X] K"^Xj6cd(Ta, n, Tc, Td), 

(Bll) 

or in Fourier space 

Xkabcd{i'^a,il^b, ifc, ifd) = 

2^ — ^ aXjbcdiit^a, it^b, IVc, IVd). (B12) 

Using this result in (B7) finally yields 

^iil^a - £a- ^aj{il^a)]X3bcd{il^a, Wb, Wc, Wd) = 
j 

+ \ ^{Uoa + UajWaf^^^iiUa, i^b, ^c, i^d)- (B13) 

3 

Subtraction of Y.jT.aj{wa)X3bcd{wa,Wb,Wc,Wd) on 
both sides of this equation and using the definition of 
Green's function, Gabi^a) ■= [{i'^a- £a)Sab- ^ab{ii^a) - 
^ab{i'^a)]~^ , we have on the left-hand-side of this equa- 
tion Y.jG~l{ii'a)X]bcd{wa,Wb,Wc,Wd)- Hencc multi- 
plying both sides of this equation by Gab{iJ^a) using ma- 
trix multiphcation and identifying the disconnected part 
of the two-particle Green's function 

Xlbcdi'i'^a, iVb, il'c, iVd) =P[Gab{iVa)Gcd{Wc)5v^,ub5v,,vi 

- Gad{Wa)Gcb{il'c)5va,i^d^iyo,iyb\^ 

(B14) 



we find for the connected part 

Xabcd{Wa, iVb, iVc, i^d) - Xl^bcdii^a, Wb, Wc, il^ d) = 

- Y Gai{il^a)^ij{i'^a)Xjbcd{i'^a, ^b, Wc, ^d) 

ij 

ij 

(B15) 

Or, using (A20), 

Xabcd{iVa,iVb,iVc,il'd) - X^abcd{i^a,il'b,iVc,il'd) = 

- ^ ^{Uij + Uji)Fl-{iUa)X3bcd{iya,iyb,iVc,i^d) 

ij 

(B16) 

Note that the inverse Green's function in (A20) cancels. 
If we had used the equation of motion for Green's func- 
tion taken with respect to its first argument, this would 
not have been the case (since in this case the inverse 
Green's function appears to the right of (A20)). 



Appendix C: Non-density-density Hamiltonians 

For the case of a general, i.e. non-density-density type 
of interaction 

Hu = Y Uijkic\c\ckCi (CI) 

ijkl 

the equation of motion involves the commutator 

[c\c]ckcu 4] = c\c]ck5ib - c\c]ci6kb- (C2) 

According to (A2), the equation of motion in this case 
generates terms which require one to measure correlation 
functions of the form 

(r.c„(r)cj(r')c](r')c,(r')). (C3) 
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